/*******************************************************************************
*
*   Copyright ©2009 Alexandru Marin
*
*   This file is part of nwave/pMRI C library.
*
*   This project started as an MSc degree project during a
*   stage at the Paris-Est university, under joint supervision of prof.
*   Bogdan Ionescu (Image Analysis and Processing Laboratory, "Politehica"
*   University of Bucharest), Caroline Chaux and Jean-Christophe Pesquet
*   (Laboratoire d'Informatique, Université Paris-Est, Institut Gaspard Monge).
*
*   nwave/pMRI C libray is free software: you can redistribute it and/or modify
*   it under the terms of the GNU General Public License as published by
*   the Free Software Foundation, either version 2 of the License, or
*   (at your option) any later version.
*
*   nwave/pMRI C libray is distributed in the hope that it will be useful,
*   but WITHOUT ANY WARRANTY; without even the implied warranty of
*   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*   GNU General Public License for more details.
*
*   You should have received a copy of the GNU General Public License
*   along with nwave/pMRI C library.  If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/


#include "mri_generator.h"

mri_image_t** mri_generate(mri_image_t* image, int reduction_factor, mri_sens_t* sens, mri_noise_cov_t* noise_cov) {

	int ncoils = sens->ncoils;
	int nslices = sens->nslices;
	int height = sens->height;
	int width = sens->width;
	int nsamples = image->nsamples;

	int newheight = height / reduction_factor;

	printf("* MRI: Generating noisy pMRI scans...\n* MRI: ncoils=%d; nsamples=%d; nslices=%d; height = %d; width = %d; newheight = %d \n\n",
			ncoils, nsamples, nslices, height, width, newheight);

	// Loop iterators
	int i, rownum, colnum, coilnum, slicenum, samplenum, roworig;

	// Allocate output images
	mri_image_t** result = (mri_image_t**) malloc(ncoils * sizeof(mri_image_t*));
	assert(result != NULL);

	for (i = 0; i < ncoils; i++) {
		result[i] = mri_image_alloc(newheight, width, nslices, nsamples);
		assert(result[i] != NULL);
	}

	// Will hold the dot product of values in the input image and corresponding sensitivity values
	complex float dotprod;

	// Will hold a vector of complex noise
	complex float* noise = (complex float*)malloc(noise_cov->size*sizeof(complex float));

	// Iterate through time dimension
	for (samplenum = 0; samplenum < nsamples; samplenum++) {

//		printf("* * * sample #%d\n", samplenum);

		// Iterate through image slices (depth)
		for (slicenum = 0; slicenum < nslices; slicenum++) {

//			printf("* * * slice #%d\n", slicenum);

			// Iterate  through the rows of the new images
			for (rownum = 0; rownum < newheight; rownum++) {

//				printf("* * * row #%d\n", rownum);

				// Iterate through columns
				for (colnum = 0; colnum < width; colnum++) {

//					printf("* * * column #%d\n", colnum);

					// NOTE: could achieve some optimization by using the same mem space
					// (i.e. not allocate noise each time we use it)
					mri_noise_cov_fill_instance(noise, noise_cov);
//					noise = mri_noise_cov_get_instance(noise_cov);

					// Iterate through the output images (each coil)
					for (coilnum = 0; coilnum < ncoils; coilnum++) {

//						printf("* * * coil #%d\n", coilnum);

						// Initialize (will accumulate a dot product)
						dotprod = 0;

						// The (rownum,colnum) pixel in the output image is given by the pixels in the
						// original image located at rows interleaved by newheight
						for (i = 0; i < reduction_factor; i++) {

							// Row number in the original image
							roworig = i * newheight + rownum;
//							roworig = (i * newheight + rownum + newheight/2 + 1) % height;

							dotprod += MRI_SENS(sens, roworig, colnum, coilnum, slicenum) * MRI_IMAGE(image, roworig, colnum, slicenum, samplenum);

//							printf("coil# %3d: (%3d, %3d) <= (%3d, %3d) \n", coilnum, rownum, colnum, roworig, colnum);
//							printf("Sens: "); mri_write_complex(stdout, MRI_SENS(sens, roworig, colnum, coilnum, slicenum));
//							printf(", Image: "); mri_write_complex(stdout, MRI_IMAGE(image, roworig, colnum, slicenum, samplenum));
//							printf(", Dotprod: "); mri_write_complex(stdout, dotprod);
//							printf("\n");

						}

//						printf("coil# %3d: (%3d, %3d) = %f %fi\n", coilnum, rownum, colnum, dotprod.re, dotprod.im);

						MRI_IMAGE(result[coilnum], rownum, colnum, slicenum, samplenum) = dotprod + noise[coilnum];

					}

				}
			}
		}
	}

	free(noise);

	return result;

}
